%% 2a
X = zeros(61*80,9);
colormap('gray');
for i = 1:9
    I = double(imread([num2str(i,'%02d') '.gif']))/256;
    subplot(3,3,i);
    X(:,i) = I(:);
    imagesc(reshape(X(:,i),61,[]));
end

m = mean(X,2);
Y = zeros(size(X,1));
for i = 1:9
    Y = Y+(X(:,i)-m)*(X(:,i)-m)';
end
Y = Y/9;

[V,D] = eig(Y);
d = diag(D);

figure;
for j = 0:3
    subplot(2,2,j+1);
    imagesc(reshape(V(:,end-j),61,[]))
    title(num2str(d(end-j)))
    axis image
    axis off
    colormap gray
end

%% 2d
I = double(imread([num2str(10,'%02d') '.gif']))/256;
J = I(:)*0;
for j = 0:3
    J = J + (I(:)'*V(:,end-j))*V(:,end-j);
end

subplot(1,2,1)
imagesc(I);
colormap gray
subplot(1,2,2)
imagesc(reshape(J,61,[]));
colormap gray
caxis([0,1])